
void MyImg::cvShiftDFT(CvArr * src_arr, CvArr * dst_arr ) {
	CvMat * tmp;
	CvMat q1stub, q2stub;
	CvMat q3stub, q4stub;
	CvMat d1stub, d2stub;
	CvMat d3stub, d4stub;
	CvMat * q1, * q2, * q3, * q4;
	CvMat * d1, * d2, * d3, * d4;

	CvSize size = cvGetSize(src_arr);
	CvSize dst_size = cvGetSize(dst_arr);
	int cx, cy;

	if(dst_size.width != size.width ||
	dst_size.height != size.height){
		cvError( CV_StsUnmatchedSizes, "cvShiftDFT", "Source and Destination arrays must have equal sizes", __FILE__, __LINE__ );
	}

	if(src_arr==dst_arr){
		tmp = cvCreateMat(size.height/2, size.width/2, cvGetElemType(src_arr));
	}

	cx = size.width/2;
	cy = size.height/2; // image center

	q1 = cvGetSubRect( src_arr, &q1stub, cvRect(0, 0, cx, cy) );
	q2 = cvGetSubRect( src_arr, &q2stub, cvRect(cx, 0, cx, cy) );
	q3 = cvGetSubRect( src_arr, &q3stub, cvRect(cx, cy, cx, cy) );
	q4 = cvGetSubRect( src_arr, &q4stub, cvRect(0, cy, cx, cy) );
	d1 = cvGetSubRect( src_arr, &d1stub, cvRect(0, 0, cx, cy) );
	d2 = cvGetSubRect( src_arr, &d2stub, cvRect(cx, 0, cx, cy) );
	d3 = cvGetSubRect( src_arr, &d3stub, cvRect(cx, cy, cx, cy) );
	d4 = cvGetSubRect( src_arr, &d4stub, cvRect(0, cy, cx, cy) );

	if(src_arr!=dst_arr){
		if( !CV_ARE_TYPES_EQ( q1, d1 )){
			cvError( CV_StsUnmatchedFormats, "cvShiftDFT", "Source and Destination arrays must have the same format", __FILE__, __LINE__ );
		}
		cvCopy(q3, d1, 0);
		cvCopy(q4, d2, 0);
		cvCopy(q1, d3, 0);
		cvCopy(q2, d4, 0);
	}
	else{
		cvCopy(q3, tmp, 0);
		cvCopy(q1, q3, 0);
		cvCopy(tmp, q1, 0);
		cvCopy(q4, tmp, 0);
		cvCopy(q2, q4, 0);
		cvCopy(tmp, q2, 0);
	}
}


int MyImg::ApplyFilter(){
	int i,j;
	double dist1,dist2,dist3,dist0;
	int fcx,fcy;

	fcx	=	(int)floor(thisWidth/4);
	fcy	=	(int)floor(thisHeight/4);


	for(i=0;i<thisWidth;i++){
		for(j=0;j<thisHeight;j++){
			dist0	=	pow(i/fcx,2)		+pow(j/fcy,2);
			dist1	=	pow((i-thisWidth)/fcx,2)+pow(j/fcy,2);
			dist2	=	pow((i-thisWidth)/fcx,2)+pow((j-thisHeight)/fcy,2);
			dist3	=	pow(i/fcx,2)		+pow((j-thisHeight)/fcy,2);
			if(dist0<1 || dist1<1 || dist2<1 || dist3<1){

			}
			else{
				setFFT(i,j,true,0);
				setFFT(i,j,false,0);
			}
		}
	}
	return(0);
}


int MyImg::ComputeFFTw(){
	IplImage	*img=NULL;
	uchar		*img_data;

	fftw_complex    *data_in;
	fftw_plan       plan_f;
	double		tempr;
	double		tempi;
	double		temp;

	double dist1,dist2,dist3,dist0;
	int fcx,fcy;



	int k, i, j;

	if(FFTexist){
	}
	else{


		//charge l'image
		img	=cvLoadImage(thisFilename.data(), CV_LOAD_IMAGE_GRAYSCALE );
		if ( img == NULL ) {
			cout << "IMAGE NOT FOUND" << endl;
			return(-1);
		}
		//note ses caractéristiques
		thisWidth  	=img->width;
		thisHeight 	=img->height;
		thisStep	=img->widthStep;

		fcx	=	(int)floor(thisWidth/4);
		fcy	=	(int)floor(thisHeight/4);
		//cout << fcx;
		//cout << fcy;

		//copie les données
		img_data	= (uchar*) img->imageData;

		data_in = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * thisWidth * thisHeight );
		fft     = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * thisWidth * thisHeight );

		//plan_f = fftw_plan_dft_1d( thisWidth * thisHeight, data_in, fft,  FFTW_FORWARD,  FFTW_ESTIMATE );
		plan_f = fftw_plan_dft_2d( thisHeight, thisWidth, data_in, fft,  FFTW_FORWARD,  FFTW_ESTIMATE );

		//copie les données dans data in
		for( i = 0, k = 0 ; i < thisHeight ; i++ ) {
			for( j = 0 ; j < thisWidth ; j++ ) {
				data_in[k][0] = ( double )img_data[i * thisStep + j];
				data_in[k][1] = 0.0;
				k++;
			}
		}
		fftw_execute( plan_f );

		ofstream fout("fft.dat");
		for( i = 0, k = 0 ; i < thisHeight ; i++ ) {
			for( j = 0 ; j < thisWidth ; j++ ) {
				//tempr				= (double)fft[k][0];
				//tempi				= (double)fft[k][1];
				tempr				= ValFFT(i,j,true);
				tempi				= ValFFT(i,j,false);
				k++;
				temp				= log10(1+pow(pow((double)tempr, 2)+pow((double)tempi, 2), 0.5));

				dist0	=	pow((double)i/fcy,2)		+pow((double)j/fcx,2);
				dist1	=	pow((double)(i-thisHeight)/fcy,2)+pow((double)j/fcx,2);
				dist2	=	pow((double)(i-thisHeight)/fcy,2)+pow((double)(j-thisWidth)/fcx,2);
				dist3	=	pow((double)i/fcy,2)		+pow((double)(j-thisWidth)/fcx,2);
				if(dist0<1 || dist1<1 || dist2<1 || dist3<1){

				}
				else{

					temp=0;
					setFFT(i,j,true,0);
					setFFT(i,j,false,0);
					//cout << ValFFT(i,j,true) << endl;


				}
				img_data[i * thisStep + j]	= (uchar)floor(temp*36.9);
				fout <<(double)temp << "\t";

			}
			fout << endl;
		}
		fout.close();

		ofstream fouto("origin.dat");
		for( i = 0, k = 0 ; i < thisHeight ; i++ ) {
			for( j = 0 ; j < thisWidth ; j++ ) {
				fouto <<(double)data_in[k][0] << "\t";
				k++;
			}
			fouto << endl;
		}
		fouto.close();
		/*
		cvNamedWindow( "original_image", CV_WINDOW_AUTOSIZE );
		cvShowImage( "original_image", img);
		cvWaitKey();
		*/
                fftw_destroy_plan(plan_f);
                fftw_free(data_in);
                //free(img_data);
                cvReleaseImage(&img);

                //fftw_free(out);

		FFTexist	=	true;
		return(0);
	}
}
double MyImg::ValFFT(int x,int y,bool IsRealPart){
    return(ValFFT(x,y,IsRealPart,GRAYSCALE_MODE));
}

double MyImg::ValFFT(int x,int y,bool IsRealPart,int mode){
	int real;
	if(y<thisWidth && x<thisHeight){
		if(IsRealPart){
			real=0;
		}
		else{
			real=1;
		}
                switch(mode){
                    case COLOR_MODE_R:
                        if(fftR==NULL){
                            return(-1);
                        }
                        return((double)fftR[x*thisWidth+y][real]);
                        break;
                    case COLOR_MODE_G:
                        if(fftG==NULL){
                            return(-1);
                        }
                        return((double)fftG[x*thisWidth+y][real]);
                        break;
                    case COLOR_MODE_B:
                        if(fftB==NULL){
                            return(-1);
                        }
                        return((double)fftB[x*thisWidth+y][real]);
                        break;
                    default:
                        if(fft==NULL){
                            return(-1);
                        }
                        return((double)fft[x*thisWidth+y][real]);
                        break;
                }
	}
}
void MyImg::setFFT(int x, int y, bool IsRealPart,double val){
    setFFT(x,y,IsRealPart,val,GRAYSCALE_MODE);
}
void MyImg::setFFT(int x, int y, bool IsRealPart,double val,int mode){
	int real;
	if(y<thisWidth && x<thisHeight){
		if(IsRealPart){
			real=0;
		}
		else{
			real=1;
		}
                switch(mode){
                    case COLOR_MODE_R:
                        fftR[x*thisWidth+y][real]	=	val;
                        break;
                    case COLOR_MODE_G:
                        fftG[x*thisWidth+y][real]	=	val;
                        break;
                    case COLOR_MODE_B:
                        fftB[x*thisWidth+y][real]	=	val;
                        break;
                    default:
                        fft[x*thisWidth+y][real]	=	val;
                        break;
                }

	}

}

int MyImg::ComputeMultiFFTw(){
    return( ComputeFFTwGeneral(COLOR_MODE_R)+
            ComputeFFTwGeneral(COLOR_MODE_G)+
            ComputeFFTwGeneral(COLOR_MODE_B)+
            ComputeFFTwGeneral(GRAYSCALE_MODE) );
}

int MyImg::ComputeFFTwGeneral(int input){
        //IplImage	*img=NULL;
	uchar		*img_data;

	fftw_complex    *data_in;
        //fftw_complex    *tempFFT;
	fftw_plan       plan_f;
	double		tempr;
	double		tempi;
	double		temp;

	double dist1,dist2,dist3,dist0;
	int fcx,fcy;



	int k, i, j;

	if(IsFFTcomputed(input)){
	}
	else{


		//charge l'image
		//img	=cvLoadImage(thisFilename.data(), CV_LOAD_IMAGE_GRAYSCALE );
                OpenImage(input);
		if ( thisimg == NULL ) {
			cout << "IMAGE NOT FOUND" << endl;
			return(-1);
		}
		//note ses caractéristiques
		thisWidth  	=thisimg->width;
		thisHeight 	=thisimg->height;
		thisStep	=thisimg->widthStep;

		fcx	=	(int)floor(thisWidth/4);
		fcy	=	(int)floor(thisHeight/4);
		//cout << fcx;
		//cout << fcy;

		//copie les données
		img_data	= (uchar*) thisimg->imageData;

		data_in = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * thisWidth * thisHeight );
                switch(input){
                    case COLOR_MODE_R:
                        fftR     = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * thisWidth * thisHeight );
                        plan_f = fftw_plan_dft_2d( thisHeight, thisWidth, data_in, fftR,  FFTW_FORWARD,  FFTW_ESTIMATE );
                        break;
                    case COLOR_MODE_G:
                        fftG     = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * thisWidth * thisHeight );
                        plan_f = fftw_plan_dft_2d( thisHeight, thisWidth, data_in, fftG,  FFTW_FORWARD,  FFTW_ESTIMATE );
                        break;
                    case COLOR_MODE_B:
                        fftB     = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * thisWidth * thisHeight );
                        plan_f = fftw_plan_dft_2d( thisHeight, thisWidth, data_in, fftB,  FFTW_FORWARD,  FFTW_ESTIMATE );
                        break;
                    default:
                        fft     = ( fftw_complex* )fftw_malloc( sizeof( fftw_complex ) * thisWidth * thisHeight );
                        plan_f = fftw_plan_dft_2d( thisHeight, thisWidth, data_in, fft,  FFTW_FORWARD,  FFTW_ESTIMATE );
                        break;
                }

		//plan_f = fftw_plan_dft_1d( thisWidth * thisHeight, data_in, fft,  FFTW_FORWARD,  FFTW_ESTIMATE );


		//copie les données dans data in
		for( i = 0, k = 0 ; i < thisHeight ; i++ ) {
			for( j = 0 ; j < thisWidth ; j++ ) {
				data_in[k][0] = ( double )img_data[i * thisStep + j];
				data_in[k][1] = 0.0;
				k++;
			}
		}
		fftw_execute( plan_f );

		ofstream fout("fft.dat");
		for( i = 0, k = 0 ; i < thisHeight ; i++ ) {
			for( j = 0 ; j < thisWidth ; j++ ) {
				//tempr				= (double)fft[k][0];
				//tempi				= (double)fft[k][1];
				tempr				= ValFFT(i,j,true,input);
				tempi				= ValFFT(i,j,false,input);
				k++;
				temp				= log10(1+pow(pow((double)tempr, 2)+pow((double)tempi, 2), 0.5));

				dist0	=	pow((double)i/fcy,2)		+pow((double)j/fcx,2);
				dist1	=	pow((double)(i-thisHeight)/fcy,2)+pow((double)j/fcx,2);
				dist2	=	pow((double)(i-thisHeight)/fcy,2)+pow((double)(j-thisWidth)/fcx,2);
				dist3	=	pow((double)i/fcy,2)		+pow((double)(j-thisWidth)/fcx,2);
				if(dist0<1 || dist1<1 || dist2<1 || dist3<1){

				}
				else{

					temp=0;
					setFFT(i,j,true,0,input);
					setFFT(i,j,false,0,input);
					//cout << ValFFT(i,j,true) << endl;


				}
				img_data[i * thisStep + j]	= (uchar)floor(temp*36.9);
				fout <<(double)temp << "\t";

			}
			fout << endl;
		}
		fout.close();

		ofstream fouto("origin.dat");
		for( i = 0, k = 0 ; i < thisHeight ; i++ ) {
			for( j = 0 ; j < thisWidth ; j++ ) {
				fouto <<(double)data_in[k][0] << "\t";
				k++;
			}
			fouto << endl;
		}
		fouto.close();
		/*
		cvNamedWindow( "original_image", CV_WINDOW_AUTOSIZE );
		cvShowImage( "original_image", img);
		cvWaitKey();
		*/
                fftw_destroy_plan(plan_f);
                fftw_free(data_in);
                //free(img_data);
                cvReleaseImage(&thisimg);
                thisimg=NULL;

                //fftw_free(out);

		FFTexist	=	true;
		return(0);
	}
}


int MyImg::ComputeFFT(){
	IplImage * im;

	IplImage * image_Im;//partie imaginare
	IplImage * image_Re;//partie reelle

	IplImage * realInput;
	IplImage * imaginaryInput;
	IplImage * complexInput;
	int dft_M, dft_N;
	CvMat* dft_A, tmp;

	double m, M;

	im = cvLoadImage( thisFilename.data(), CV_LOAD_IMAGE_GRAYSCALE );//charge l'image
	if( !im )
		return -1;

	realInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);//partie reelle
	imaginaryInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 1);//partie imaginaire
	complexInput = cvCreateImage( cvGetSize(im), IPL_DEPTH_64F, 2);//complexe

	cvScale(im, realInput, 1.0, 0.0);//met la meme echelle au reel qu a l' image?
	cvZero(imaginaryInput);//zero dans l'imaginaire
	cvMerge(realInput, imaginaryInput, NULL, NULL, complexInput);//copie les chan. non null de  realinput,imaginary dans complex

	dft_M = cvGetOptimalDFTSize( im->height - 1 );//renvoie la taille optimale pour calculer la FFT
	dft_N = cvGetOptimalDFTSize( im->width - 1 );

	dft_A = cvCreateMat( dft_M, dft_N, CV_64FC2 );//cree la matrice M*N
	image_Re = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);//cree image relle
	image_Im = cvCreateImage( cvSize(dft_N, dft_M), IPL_DEPTH_64F, 1);//cree image imaginaire

	// copy A to dft_A and pad dft_A with zeros
	cvGetSubRect( dft_A, &tmp, cvRect(0, 0, im->width, im->height));//comme taille!=taille optimale, il coupe/copie
	cvCopy( complexInput, &tmp, NULL );
	if( dft_A->cols > im->width ) {
		cvGetSubRect( dft_A, &tmp, cvRect(im->width, 0, dft_A->cols - im->width, im->height));
		cvZero( &tmp );
	}

	// no need to pad bottom part of dft_A with zeros because of
	// use nonzero_rows parameter in cvDFT() call below

	cvDFT( dft_A, dft_A, CV_DXT_FORWARD, complexInput->height );

	cvNamedWindow("win", 0);
	cvNamedWindow("magnitude", 0);
	cvShowImage("win", im);

	// Split Fourier in real and imaginary parts
	cvSplit( dft_A, image_Re, image_Im, 0, 0 );

	//******************
	cvPow( image_Re, image_Re, 2.0);
	cvPow( image_Im, image_Im, 2.0);
	cvAdd( image_Re, image_Im, image_Re, NULL);
	cvPow( image_Re, image_Re, 0.5 );

	// Compute log(1 + Mag)
	cvAddS( image_Re, cvScalarAll(1.0), image_Re, NULL ); // 1 + Mag
	cvLog( image_Re, image_Re ); // log(1 + Mag)


	// Rearrange the quadrants of Fourier image so that the origin is at
	// the image center
	cvShiftDFT( image_Re, image_Re );

	cvMinMaxLoc(image_Re, &m, &M, NULL, NULL, NULL);
	cvScale(image_Re, image_Re, 1.0/(M-m), 1.0*(-m)/(M-m));
	cvShowImage("magnitude", image_Re);
	//****************************

	//****************************
	cvWaitKey();
	//cout << dft_A->data.ptr << endl;
	FFTexist	=	true;
	cvReleaseImage(&im);
	cvReleaseImage(&image_Im);
	cvReleaseImage(&image_Re);
	cvReleaseImage(&realInput);
	cvReleaseImage(&imaginaryInput);
	cvReleaseImage(&complexInput);
	return(0);

	/*IplImage * im;

	IplImage * image_Im;//partie imaginare
	IplImage * image_Re;//partie reelle

	IplImage * realInput;
	IplImage * imaginaryInput;
	IplImage * complexInput;
	 **/
}


bool MyImg::IsFFTcomputed(int mode){
    switch(mode){
        case COLOR_MODE_R:
            if(fftR==NULL){
                return(false);
            }
            else{
                return(true);
            }
        case COLOR_MODE_G:
            if(fftG==NULL){
                return(false);
            }
            else{
                return(true);
            }
        case COLOR_MODE_B:
            if(fftB==NULL){
                return(false);
            }
            else{
                return(true);
            }
        default:
            if(fft==NULL){
                return(false);
            }
            else{
                return(true);
            }
    }
}

